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We present a universal technique for quantum state esti- 
mation based on the maximum-likelihood method. This ap- 
proach provides a positive definite estimate for the density 
matrix from a sequence of measurements performed on iden- 
tically prepared copies of the system. The method is versatile 
and can be applied to multimode radiation fields as well as 
to spin systems. The incorporation of physical constraints, 
which is natural in the maximum-likelihood strategy, leads to 
a substantial reduction of statistical errors. Numerical imple- 
mentation of the method is based on a particular form of the 
Gauss decomposition for positive definite Hermitian matrices. 

PACS Numbers: 03.67, 03.65.Bz 

In quantum mechanics, the achievable information on 
a physical system is encoded into the density matrix g, 
which allows one to evaluate all possible expectation val- 
ues through the Born statistical rule (O) — Tr (gO). In 
order to obtain full information on a quantum system we 
need to estimate its density matrix. In principle, this 
can be accomplished by successive measurements on re- 
peated identical preparations of the same system. With 
a proper choice of the measurements, and after collecting 
a suitably large number of data, we can arrive at reliable 
knowledge of the quantum state of the system. 

The problem of inferring the complete quantum state 
from experimental data has received a lot of attention 
over past several years. Physical systems whose quan- 
tum state has been fully characterized in recent exper- 
iments, include now a single light mode [jjj, a diatomic 
molecule j2j, a trapped ion ||, and an atomic beam jjj. 
These fascinating advances stimulate further theoretical 
research in two main directions: on one hand, in imple- 
menting effective measurement schemes that connect the 
density matrix to directly observable quantities. On the 
other hand, in designing efficient data processing algo- 
rithms in a practical experimental setup in order to ex- 
tract the optimal amount of information on the quantum 
state. In a laboratory, we always deal with finite ensem- 
bles of copies of the measured system || . In addition, 
the process of detection is usually affected by various im- 
perfections. This implies the need of developing novel 
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tools specifically designed to process realistic and finite 
experimental samples. 

In this Communication we present a general-purpose 
method for quantum state estimation based on the 
maximum-likelihood (ML) approach ||. We consider 
statistical treatment of a sample of measurements per- 
formed on repeated preparations of a given system. The 
approach presented in this Communication is very gen- 
eral: it allows one to extract the information on the quan- 
tum state from data collected in a generic scheme, with- 
out assuming any specific form of the measurement. Its 
principle of operation is to find the quantum state that 
is most likely to generate the observed data. This idea 
is quantified and implemented using the concept of the 
likelihood functional. 

The ML strategy is an entirely different approach to 
quantum state measurement compared to the standard 
quantum-tomographic techniques (?]||]. In quantum to- 
mography the expectation value of an operator is ob- 
tained by averaging a special function (so called "pat- 
tern function" ) of experimental data of a sufficiently com- 
plete set of observables — a "quorum" of observables. In 
homodyne tomography the quorum observables are the 
quadratures of the e.m. field for varying phase with re- 
spect to the local oscillator. Hence, typically, a matrix 
element of the quantum state is obtained by averaging 
its pertaining pattern function over data. This method 
is very general and efficient, however, in the averaging 
procedure, the matrix elements are allowed to fluctuate 
statistically through negative values, with resulting large 
statistical errors. 

In contrast, the ML method estimates the quantum 
state as a whole. Such a procedure incorporates a priori 
knowledge about relations between elements of the den- 
sity matrix. This guarantees positivity and normaliza- 
tion of matrix, with the result of a substantial reduction 
of statistical errors. These advantages of the ML ap- 
proach are inevitably related to increased computational 
complexity of the estimation procedure, which remains a 
highly nontrivial problem even if we resort to numerical 
means. To the best of our knowledge, we present in this 
Communication the first general solution to this prob- 
lem, which provides an effective numerical algorithm for 
the ML estimation of the density matrix. 

We start with the derivation of the likelihood func- 
tional £(g), which links the raw experimental results with 
the object to be reconstructed, i.e. the density matrix. 
The physical situation we have in mind is an experi- 
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ment consisting of iV measurements performed on iden- 
tically prepared copies of a given system. Quantum me- 
chanically, each measurement is described by a positive 
operator- valued measure (POVM). The outcome of the 
ith measurement corresponds to the realization of a spe- 
cific element of the POVM used in the corresponding run. 
We shall denote this element by Ti . The likelihood func- 
tional £(§) describes the probability of obtaining the set 
of outcomes for a given density matrix g. For measure- 
ments performed on repeated preparations of the system, 
it is given by the product 



(i) 



After the experiment is performed, the operators T% are 
determined by the outcomes of the measurements. The 
unknown element of the above expression, which we want 
to infer from our data, is the density matrix describing 
the measured ensemble. The general estimation strat- 
egy of the ML technique is to maximize the likelihood 
functional over the set of the density matrices. Several 
properties of the likelihood functional are easily found, if 
we restrict ourselves to finite dimensional Hilbert spaces. 
In this case, it can be easily proved that C(g) is a con- 
cave function defined on a convex and closed set of den- 
sity matrices. Therefore, its maximum is achieved either 
on a single isolated point, or on a convex subset of den- 
sity matrices. In the latter case, the experimental data 
are insufficient to provide a unique estimate for the den- 
sity matrix using the ML strategy. On the other hand, 
existence of a single maximum allows us to assign unam- 
biguously the ML estimate for the density matrix. This 
estimate satisfies all the physical constraints, such as nor- 
malization and positivity. 

ML estimation of the quantum state, despite its el- 
egant general formulation, presents a highly nontrivial 
constrained optimization problem, even if we resort to 
purely numerical means. The central difficulty lies in 
the appropriate parameterization of the set of all density 
matrix. The parameter space should be of the minimum 
dimension in order to preserve the maximum of the like- 
lihood function as a single isolated point. Additionally, 
the expression of quantum expectation values in terms 
of this parameterization should enable fast evaluation of 
the likelihood function, as this step is performed many 
times in the course of numerical maximization. 

Here, we introduce a parameterization of the set of 
density matrices which provides an efficient algorithm for 
maximization of the likelihood function. We represent 
the density matrix in the form 
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(2) 



which automatically guarantees that g is positive and 
Hermitian. The remaining condition of unit trace Tig = 
1 will be taken into account using the method of La- 
grange multipliers. In order to achieve the minimal pa- 



rameterization, we assume that T is a complex lower tri- 
angular matrix, with real elements on the diagonal. This 
form of T is motivated by the Cholesky decomposition 
known in numerical analysis J^] for arbitrary non neg- 
ative Hermitian matrix. For an Af-dimensional Hilbert 
space, the number of real parameters in the matrix T 
is M + 2M(M - l)/2 = A/ 2 , which equals the number 
of independent real parameters for a Hermitian matrix. 
This confirms that our parameterization is minimal, up 
to the unit trace condition. 

In numerical calculations, it is convenient to replace 
the likelihood functional by its natural logarithm, which 
of course docs not change the location of the maximum. 
Thus the function subjected to numerical maximization 
is given by 
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L(f) = ^lnTrfTtJTi) - ATr(f' t f) 
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where A is a Lagrange multiplier accounting for normal- 
ization of g that equals the total number of measurements 
N jnj . This formulation of the maximization problem 
allows one to apply standard numerical procedures for 
searching the maximum over the M 2 real parameters of 
the matrix T. The examples presented below use the 
downhill simplex method |12|. 

Our first example is the application of the ML estima- 
tion in quantum homodyne tomography of a single-mode 
radiation field Q], which is so far the most successful 
method in measuring nonclassical states of light Ji],pT|. 
The experimental apparatus used in this technique is the 
homodyne detector. The realistic, imperfect homodyne 
measurement is described by the positive operator- valued 
measure 
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where r\ is the detector efficiency, and x v is the quadra- 
ture operator, depending on the externally adjustable lo- 
cal oscillator (LO) phase <p. 

After repeating the measurement N times, we obtain 
a set of pairs (xi;tpi) consisting of the outcome Xi and 
the LO phase ip% for the zth run, where i = 1,. . . , JV. 
The log-likelihood functional is given by Eq. (||) with 
T = H(xi-(fi). Of course, for a light mode it is neces- 
sary to truncate the Hilbert space to a finite dimensional 
basis. We shall assume that the highest Fock state has 
M — 1 photons, i.e. that the dimension of the truncated 
Hilbert space is M. For the expectation Tr[T^TH{x; ip)} 
it is necessary to use an expression which is explicitly 
positive, in order to protect the algorithm against oc- 
currence of small negative numerical arguments of the 
logarithm function. A simple derivation yields 
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H n (x) exp(— x 2 /2)/V2 n n\ir 1 / 2 are eigenstates of the har- 
monic oscillator in the position representation — H n (x) 
being the nth Hermite polynomial. 

We have applied the ML technique to reconstruct the 
density matrix in the Fock basis from Monte Carlo simu- 
lated homodyne statistics. Fig. [j] depicts the matrix ele- 
ments of the density operator as obtained for a coherent 
state and a squeezed vacuum, respectively. Remarkably, 
only 50000 homodyne data have been used for quantum 
efficiency at photodetectors ra = 80%. 

Since statistical aspects of standard quantum homo- 
dyne tomography have been thoroughly studied [p"3| , this 
gives us an opportunity to compare it with the ML esti- 
mation. In the tomographic approach, statistical errors 
are known to grow rapidly with decreasing efficiency ra 
of the detector. In contrast, the elements of the den- 
sity matrix reconstructed using the ML approach remain 
bounded, as the whole matrix must satisfy positivity and 
normalization constraints. This results in much smaller 
statistical errors. As a comparison one could see that 
the same precision of the reconstructions in Fig. [I] could 
be achieved using 10 7 -10 8 data samples with the con- 
ventional quantum tomography of Ref. Jt|] . On the other 
hand, in order to find numerically the ML estimate we 
need to set a priori the cut-off parameter for the photon 
number, and its value is limited by increasing computa- 
tion time. 

Another relevant example is the reconstruction of the 
quantum state of two-mode field using single-LO homo- 
dyning [ p"4[ . Here, the full joint density matrix can be 
measured by scanning the quadratures of all possible lin- 
ear combinations of modes. For two modes the measured 
quadrature operator is given by Xe^ $ x — (ae~ l ^° cos 9 + 
be'^ 1 sin 9 + h.c.)/V2, where (0,V>o,^i) € S 2 x [0,2tt], 
S 2 being the Poincare sphere and one phase ranging be- 
tween and 2ir. In each run these parameters are chosen 
randomly. The POVM describing the measurement is 
given by the right-hand side of Eq. (^) , with x v replaced 
by &0i/>oV>i' ano - the quantum expectation values of the 
POVM can be written as 
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We have simulated an experiment for the two orthogo- 
nal states = (| 00) + |11))/V2 and |* 2 ) = (|01) + 
|10))/-\/2- We reconstructed the density matrix in the 
two-mode Fock basis using the ML technique. The re- 
sults are depicted in Fig. g. 

Finally, we mention that the ML procedure can be ap- 
plied also for reconstructing the density matrix of spin 



systems. For example, let us consider N repeated prepa- 
rations of a pair of spin-1/2 particles. The particles are 
shared by two parties. In each run, the parties select 
randomly and independently from each other a direction 
along which they perform spin measurement. The ob- 
tained result is described by the joint projection operator 
(spin coherent states) 



(J) 



where ttf and Of are the vectors on the Bloch sphere 
corresponding to the outcomes of the ith run, and the 
indices A and B refer to the two particles. As in the 
previous examples, it is convenient to use an expression 
for the quantum expectation value Tr(T>TjFi) which is 
explicitly positive. The suitable form is 



Tr(ftf^)=^|(/.|f|Qf,Of)| 2 
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where is an orthonormal basis in the Hilbert space of 
the two particles. The result of a simulated experiment 
with only 500 data for the reconstruction of the density 
matrix of the singlet state is shown in Fig. g. 

We conclude this Communication with a brief discus- 
sion of the statistical uncertainty of the ML estimate. 
The likelihood function can be formally regarded as a 
probability distribution on the parameter space. In our 
case, this space is spanned by M 2 real parameters which 
constitute the triangular matrix T. We shall denote these 
parameters in the vector form as t. The formal distri- 
bution is given, up to the normalization constant, by 
<5[Tr(T tf ) - 1] exp L(f). In the limit of the large number 
of measurements, expi(T) takes the form of the Gaus- 
sian Jljj ], with the quadratic form in the exponent given 
by the matrix G — —d 2 L/dtdt'. Furthermore, the con- 
straint Tr(TtT) = 1 means locally orthogonality to the 
gradient u = <9Tr(T^T)/<9t. The covariance matrix for 
the parameters t is consequently given by |i~6[ ] 



V = G~ 



U T G % 
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With this result, we can estimate errors for the den- 
sity matrix using simply the propagation law applied to 
Eq. (g). 

Summarizing, we have developed a universal maxi- 
mum likelihood algorithm for estimating the density ma- 
trix. With respect to conventional quantum tomography 
this method has the great advantage of needing much 
smaller experimental samples, making experiments with 
low data rates now feasible, however with a truncation 
of the Hilbert space dimension. We have shown that the 
method is general and the algorithm has solid method- 
ological background, its reliability being confirmed in a 
number of Monte Carlo simulations. 
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FIG. 1. ML reconstruction of the density matrix of a sin- 
gle-mode radiation field. On the left the matrix elements 
obtained for a coherent state with {a^a) = 1 photon. On the 
right for a squeezed vacuum with (a^a) — 0.5 photon. In 
both cases the ML technique has been applied to a sample of 
50000 simulated homodyne data, and for quantum efficiency 
77 = 80%. 





FIG. 2. ML reconstruction of the density matrix of a 
two-mode radiation field. On the left the matrix elements 
obtained for the state = (jOO) + |ll))/\/2; on the right 
for |*2) = (|01) + |10})/\/2. For we used 100000 simu- 
lated homodyne data and n — 80%; for we used 20000 
data and n = 90%. 




FIG. 3. ML reconstruction of the density matrix of a pair of 
spin- 1/2 particles in the singlet state. The particles are shared 
by two parties. In each run, the parties select randomly and 
independently from each other a direction along which they 
perform spin measurement. The matrix elements has been 
obtained by a sample of 500 simulated data. 
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